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Abstract 

Under some weak conditions, the first-passage time of the Brownian 
motion to a continuous curved boundary is an almost surely finite stopping 
time. Its probability density function (pdf) is explicitly known only in few 
particular cases. Several mathematical studies proposed to approximate 
the pdf in a quite general framework or even to simulate this hitting time 
using a discrete time approximation of the Brownian motion. The authors 
study a new algorithm which permits to simulate the first-passage time 
using an iterating procedure. The convergence rate presented in this paper 
suggests that the method is very efficient. 

Key words and phrases: first-passage time, Brownian motion, potential the¬ 
ory, randomized algorithm. 
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Introduction 

Modeling biological or physical systems often requires handling one-dimensional 
diffusion processes. The marginal probability distribution of such processes, at 
a fixed time, permits a quite precise description of the model. Nevertheless, 
in many applications, this information is insufficient and the description of the 
whole path becomes crucial. This is namely the case for a variety of problems 
related to neuronal sciences, financial derivatives with barriers, ruin probability 
of an insurance fund, optimal stopping problems,... In these frameworks, the 
main task is the description of the first passage time densities for time-dependent 
boundaries. Let us just mention some references in engineering reliability 
epidemiology [40j, biology [34], mathematical finance [18, (36) 33 and references 
concerning the framework of level-crossing problems mm- 

For instance, let us focus our attention on a simple interpretation of neural 
transmission. When a neuron is stimulated by pressure, heat, light, or chemical 
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information, its membrane voltage changes as time elapses and, as soon as it 
reaches a constant threshold, the depolarization phenomenon occurs and the 
voltage is reset to a resting potential. The family of integrate-and-fire spiking 
neuron models is based on this simple interpretation. The firing time therefore 
corresponds to the first-passage time of the membrane potential, represented by 
a stochastic mean-reverting process (usually the Ornstein-Uhlenbeck process) to 
the neural threshold (Giorno et al. [20], Lansky et al. [25], Wan and Tuckwell 
m , for an introduction to noise in the nervous system see Part I Chapter 5 in 
[15] , for the integrate-and-fire model see Chapter 10 in [15]). 

Our main motivation is to emphasize an algorithmic approach in order to 
approximate the first-passage time of the Brownian motion to curved bound¬ 
aries. The field of application of such an algorithm at a first glance may appear 
as quite restrictive since it concerns the Brownian motion but in fact a lot of 
families of diffusion processes are concerned. Indeed it is possible to express 
various stochastic paths as functions of Brownian paths in the spirit of Wang 
and Potzelberger [42[. Hence using simple time transformations, we are going 
to present an application of the results to the Ornstein-Uhlenbeck process (see 
Section [3]) . 

In order to describe approximations of the first-passage time of the Brownian 
motion, we assume that this stopping time is almost surely finite. In this way, we 
introduce particular conditions for this property to be satisfied. Let us consider 
a continuous function ip : R + —> R. satisfying the following hypothesis: 

<p(0) > 0 and limsup — < 1. (HI) 

t-i- oo v2tloglogt 

We then define the hitting time 

= inf{t > 0 : B t = ip(t)} (0.1) 

where (H t , t > 0) stands for a standard one-dimensional Brownian motion. 
Under ( |H1| ), the a.s. finiteness of t v is an obvious consequence of the law of the 
iterated logarithm (see e.g.[28] Th.9.23 p.112]). It is quite difficult to obtain 
precise information about this stopping time in general situations. 

The study of the approximation of the hitting times for Brownian motion 
and general Gaussian Markov processes is an active area of research. Several 
alternatives for dealing with the characterization of hitting times exist. 

Approximation of the probability density function 

For particular cases, the probability density function of the Brownian passage 
time can be computed explicitly. Lerche m used the method of images in order 
to obtain explicit expressions of the p.d.f. p defined by p{t) dt = P (t v £ dt ). 
However only few cases are concerned by such a study. 

Durbin puna proposed to approximate the first-passage distribution p(t) 
of the Brownian motion as follows: p can be represented by an expansion 

k 

P(t) = + (—1 ) k r k {t), k> 1, 

1=1 

where qj for 1 < j < k and r}. are defined by multiple integrals depending on the 
boundary <p. The approximation simply consists in truncating the expansion. 
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Let us note that the first term corresponds in fact to the tangent approximation 
of Strassen [3S] and Daniels [5] , see also m- The convergence of the series and 
the error bounds can be made precise if the curved boundary is wholly concave 
or wholly convex. Many studies concern such a series expansion, let us mention 
a few: Ferebee El, Ricciardi et al. [35]; Giorno et al. EH, Sacerdote and 
Tomassetti E to deal with more general diffusion processes. The numerical 
approach proposed in JB] seems to be particularly efficient. In E, the authors 
proposed a comparison between the approximation developed by Durbin and an 
other numerical resolution of the Volterra equation for Gaussian processes. 

One method in approximating the passage time of a Brownian motion, or 
even of a quite general diffusion, through a curved boundary is to replace the 
initial boundary by an other one which is close and which leads to an explicit 
expression of the hitting time probability. Such method permits to obtain some 
bounds. It was first introduced for the Brownian motion in [3], and applied 
for instance to piecewise continuous boundaries [42]. Finally an other method 
consists in writing the p.d.f of the hitting time as the expectation of a par¬ 
ticular functional of a three-dimensional Brownian bridge. It suffices then to 
approximate this expectation through a Monte Carlo method ]2Bj. 

Approximation of the first passage time. 

All methods described so far concern the approximation of the pdf. It can be of 
particular interest to simulate directly the first passage time t v or to compute 
the probability for the hitting time to be smaller than some given T > 0, without 
computing the pdf. The solution consists in using a time discretization of the 
Brownian motion on [0, T). The time interval is then split into n small intervals 
of the kind [(k — 1 )T/n,kT/n\, with 1 < k < n. It is therefore possible just 
to simulate the hitting time of the corresponding Euler scheme. Of course this 
should upper-bound the stopping time. One solution to overcome this problem 
is to improve the algorithm by shifting the boundary to reach: we stop the Euler 
scheme as soon as it exits from a suitable smaller domain. Let us note that this 
general procedure can also be applied to diffusion processes. It has been first 
introduced for geometrical Brownian motion (finance) in [5] and then extended 
to general diffusions with nice coefficients in [251. 

An other method in order to improve the approximation of the hitting time 
consists in testing, at each endpoint kT/n, if the event B^T/n < <p(kT/n) is sat¬ 
isfied and if the Brownian path on the small intervals, conditionally on its value 
at the end point, hits the curved boundary. This method can also be applied 
to diffusions and needs therefore precise asymptotics of hitting probabilities for 
pinned diffusions. A first important study in that direction is [23] where the 
coefficients of the diffusion are frozen at the starting point on each small in¬ 
terval leading to asymptotics of the probabilities. Nevertheless the method can 
become onerous if the observed time interval [0, T\ is large and sometimes gives 
incorrect asymptotics: it has been pointed out, by the numerical treatment of 
some precise examples, that the approximations produced by this method can 
be far from the true ones. See for this point Giraudo and Sacerdote [22] (O.U. 
process and Feller model), who also suggest some formulas for the computation 
of the crossing probability, see also [23] . Baldi and Caramellino [2] presents 
precise asymptotics for general pinned diffusions which permits to improve the 
approximation of hitting times. Such results can be developed further in the 
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particular gaussian framework for any dimension [8]. 

A new algorithm. 

The aim of this study is to present a new method of approximation of t v . Let 
us explain intuitively the simulation procedure. If ip is an increasing curve with 
cp(0) > 0, then the Brownian motion needs to successively cross a sequence of 
imaginary horizontal lines before hitting the boundary. The first line to cross 
corresponds to the value </?(0) and needs a random time denoted by 7i- At that 
time, the value of the curved boundary is ip{ 71). The Brownian motion therefore 
needs to cross this second horizontal line, it shall happen at time I 2 and the 
new horizontal line to cross becomes 99 ( 72 ) and so on... Figure [l] (left) illustrates 
this procedure. The sequence of stopping times (' T n ) converges towards t v and 
will be used in order to obtain an approximation. We shall introduce a stopping 
procedure in this sequence of random times which depends on a small parameter 
e associated to the error size of the approximation: the sequence is stopped as 
soon as the distance between two successive horizontal lines is smaller than 
e. The outcome of the algorithm corresponds therefore to a random variable 

which can be exactly simulated and such that r<£ converges toward t v in 
distribution as e tends to 0. 

The algorithm can be modified when the curved boundary ip does not satisfy 
the monotonic property anymore. In such a slightly different context, it suffices 
to tilt the successive imaginary horizontal lines in such a way that the common 
slope corresponds to inf t >o see Figure [l] (right). 



Figure 1: Illustration of the algorithm for an increasing boundary ip with its 
associated successive horizontal lines (left) and for a general boundary (right). 

To sum up, two different families of sequences will be developed and the 
associated convergence rates are estimated. The first algorithm developed in 
Sectionconcerns increasing curved boundaries and the second one, Section [2] 
permits us to deal with quite general boundaries provided that its derivative is 
bounded. In the last section, we present different examples in order to illustrate 
the algorithm efficiency. 
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1 First-passage time to non-decreasing bound¬ 
aries 


Let us assume that the boundary ip satisfies (HI) and that the following addi¬ 
tional conditions hold 


ip : R + — > ffi. is a non-decreasing (^-continuous function, 

2 <p'(t)VT+t < 1 , \/t> 0 . 


(H2) 

(H3) 


We introduce the algorithm associated to the hitting time t v defined by (0.1). 


Algorithm (Al). Let e > 0 be a small parameter and (G n ) n >o fl sequence of 
independent standard Gaussian distributed random variables. 

Initialization: T 0 = 0, T\ = (<p(0)/Go) 2 and Af t = 1. 

While p{Ti) — p(Tq) > e do: 


(T 0 , TO <- (Ti,Ti + - ^(T 0 )) 2 /G^) 

Af e <-AT t + 1 . 


Outcome: •<— T\ and Af e . 

Let us just note that Algorithm |(A1)| is very simple to use since each step 
only requires one Gaussian distributed random variable. Moreover it is a ap¬ 
proximation of the first-passage time: 


Theorem 1.1. 

(H2| and 


1. Let us assume that the boundary function satisfies EL 


(H3| 


then the random variable rf n defined in Algorithm (Al) 


converges in distribution towards t v 
M ore precisely 


defined by (0.11 as e tends to zero. 


F e (t - e) - ^ = < F(t) < F e (t), for any t > e, (1.2) 

V Z7T 

where F (resp. F e ) is the cumulative distribution function of t v (resp. 

T V- 

2. There exists a constant C > 0 such that the random number of iterations 
Af e defined in Algorithm \ (A 1)\ satisfies: 

E[A4] < G^ObgO- (1-3) 


The parameter e describes the precision of the approximation. The number 
of steps in the Algorithm |(A1)| is very small (even smaller than usual results 
obtained for algorithms based on random walks on spheres, which are close 
to Algorithm |(A1)| see [32]) : in fact the constant appearing in can be 
explicitly computed: for any constant 0 < k < 1/2, there exists Cq(k) > 0 such 


that (1.3) is satisfied as soon as e < eg, with the particular constant 


G=—, m = log(4) + ^3 p 

rriK a/'k 


and 


poo 

/i = / (log\x\)e~ x dx. 

Jo 


(1.4) 
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The proof of Theorem |l.l| is based on a main argument developed in the following 
proposition: each step of Algorithm |(A1)| has to be related to a particular part 
of the Brownian paths before hitting the boundary. 

Proposition 1.2. Let ( B t , t > 0) be a standard one-dimensional Brownian 
motion. We define the following sequence of stopping times: So = To = 0 and 
for any n > 1; 


:= inf 


|t>0: Bt+Tn-x = ¥>(7^-i)j and T n ■= si + • ■ • + s n , (1.5) 


where the function ip satisfies (HI I, (H2l and (H3). Then the following proper¬ 
ties hold: 


1. (Tn)n>o is a non-decreasing sequence which almost surely converges to¬ 
wards T,p. 

2. Let n > 1, then the probability distribution of s n+ 1 given the a-algebra 
T n := tr{si...,s„} is identical as ( <p(T n ) - tp(T n -i)) 2 /G 2 where (G„)„> 0 
is a sequence of independent standard Gaussian random variables. More¬ 
over si = (tp(0)/Go) 2 . 


3. Let A4 e := inf{n > 1 : ip(l~ n ) — Bj- n < e}, then and t^, defined in 
Algorithm\(Al j\ are identically distributed, so are and Af e . 


Let us note that the mean of each random variable s n defined by (1.51 is infi¬ 
nite since E[G -2 ] = +00 where G is a standard Gaussian variable. Proposition 


1.2 suggests that the first-passage time can be obtained as a sum of positive 


random variables of infinite average, we easily deduce E[r ¥ ,] = + 00 . In the 
particular case of increasing boundaries <p, the sum has infinitely many terms. 


Proof of Proposition \1.2\ 

Step 1. By construction, the sequence (7^) n >o is non-decreasing and non¬ 
negative: it converges almost surely to Too. Since ip is a non-decreasing bound¬ 
ary, T n < t v for any n > 0. In particular Tx> is less than t v which is a finite 
stopping time due to the law of the iterated logarithm, see (HI I followed by 
discussion. Consequently, the random variable Bj-^ is well defined. Since p> is 
non-decreasing, we get Bj- n = ip(T n -i) for any n > 1. Taking the large n limit 
leads to Bj^ = plfToo ), the Brownian paths and the function p being continu¬ 
ous. We deduce that Too = Tp. 

Step 2. Let us first consider the stopping time Si. Using the reflection principle 
of the Brownian paths and a scaling property, we obtain: 


P(si > t) = PI sup B u < 

' 0<u<t 


=n\B t \<pm 


= P (B 2 < tp(0) 2 /t) = P(<p(0) 2 /Gq > t), t> 0. 


The general n-th case can be proven using similar arguments combined with the 
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Markov property of the Brownian motion: 

P(s n+ i > t\T n ) = P( sup B u < ip(T n ) 

y Tn.<U<T n +t 

= p( sup B u+Tri - B Tn < < 

v 0 <u<t 

= p( sup B u < <p{Tn) - tp( 

v 0<u<t 

where B is a Brownian motion independent of F n . 

Step 3. Using the results developed in Step 2, we observe that (s n ) nA ^ 4 E and 
the sequence of values Ti, defined in Algorithm |(Al)| have the same distribution. 
It is therefore obvious that Tm, an d r* are identically distributed. Indeed the 
stopping time can be rewritten as follows: 



M e = inf{?r > 1 : ip(T n ) - p(T n -i ) < e}. 


Proof of Theorem 

Step 1. Let us recall that T n is defined by (1.5). By Proposition 
for any n > 0 and in particular 7 ~m c < T ip- Hence 


1.2 


( 1 . 6 ) 

□ 


T < r 

> n — 1 ip 


P(‘ Tm t <t)> < t), Vt > 0. 

Since has the same distribution as Tm we obtain 

F&) > F(t), Mt > 0, (1.7) 

where F e and F are the associated cumulative distribution functions. Let us 
now prove the second bound in (1.2 1. For t > e, 

F e {t - e) = P(r* <t-e)= P(7m 6 <t-e) 

< P(Tm 5 < t - e, > t) + P (t v < t) 

< P(|TM e - T v \ > e) + F(t). (1.8) 

Combining the Markov property of the Brownian motion and the reflection 
principle leads to 

P e := P(|Tm c - t v \ > e) < 1 - P( sup B TMe+u > sup </?(Tm« + uj) 

<u<e 0<u<e ' 

”( sup B TMe+u -B TMe > sup <^(7Ai e + u) - ip{T Me ) + e) 

' 0<n<e 0<n<e ' 

"( sup B r Me +u - B t > <p{T Mt + e) - P(T M J + e) 


< 1 
< 1 


< 1- 


■( 


0<1fc<€ 

sup B u > 


' 0<n<€ 

< 1-PMBJ > 


sup <p'(0)e + e 
7"]\4 e “l _e 


sup ip'(6)e + e 

l~A4 e “I -6 


)■ 


7 










Using Hypothesis (H3l and straightforward computations permits us to obtain 


P(| T m . -t v | > e) < 1 - P(|B e | > 3e/2) < 3 


(1.9) 


The lower bound in (1.21 holds due to both (1.8) and (1.9). 

Step 2. Let us now focus our attention to the efficiency of this algorithm. We 
need to estimate the number of steps which depends on the small parameter e. 
Using the third result presented in Proposition 1.2 on one hand and (1.6 1 on 
the other hand, we obtain 

F(AT e > n) = F(M e > n) = P(<p(7i) - ^(7o) <p(T n ) - <p(T n -i) > e)- 


Hypothesis (H31 implies 

P(A4 > n) < P(si >2 e,... ,s n > 2e ). 


( 1 . 10 ) 


Step 2.1. Let us first estimate the previous upper-bound. We introduce a se¬ 
quence of independent standard Gaussian random variables {G n ) n >o and define 


X n = log(4G 2 ), 


= ^2 Xk and Zn = ^2 ‘ 


(i.ii) 


k =0 


k =0 


Let us define II(n, e) := P(s n > 2e). By Proposition 1.2, we know that the 
random variables 5 n +i are related to G n and therefore 

n(l,e) =P(2eG 2 <^(0) 2 ) =p(log(4G 2 ) < -log(e) + log(2) + 2log 

= f(z 0 < — log(e) + log(2) + 2 log ^(0)^ . 

Let us prove that, for n > 1, we have the general formula: 

n(n, e) < P(^„_i < - log(e) + (2 n - 1) log(2) + (2 n) log 

By Proposition 1 1.2 [ we have for n > 2, 


( 1 . 12 ) 


H(n, e) = p((^(T„_i) - <p(r„_ 2 )) 2 > 2eG 2 ^ 


(1.13) 


Since <p is a non decreasing function satisfying Hypothesis (H3), the following 
upper-bound holds for n > 2: 


<p(T n -i) - <p(T n -*) < r "~ 1 Tn 2 < 


Sn—1 


2y/l + Tn—2 2-y/l + S n —2 


(1.14) 


Hence for n = 2, (1.13) and (1.14) imply 

„2 


H(2, e) < p(|| > 2eG 2 ) = p(e(2G?)(2G 2 ) 2 < ^(0) 4 

= P(2X 0 + X-! < - log(e) + 3 log(2) + 4 log <p(0)) 
= P {Zi < -log(e) + 31og(2)+41og^(0)). 

















Using the lower-bound l+s ra _i > s ra _i and similar arguments as those developed 
previously, the general case is expressed as follows: 


II(n, e) < 


’ra-l 


2 2 (1 +«„- 2 ) >2eG ”-' " 




D n—3 


n n— 1 


< 


< 


< 


> e22 2 2 4 ... 2 2 I"- 2 >g;_,G*_ 2 ... G 2 *- 2 ') 


•(( 4 ?: 


1-1 1 


> 


e22 2 2 4 ... 2 2 ("- 2 )G 2 _ 1 G 4 _ 2 ... G^- 2) ) 


■ > e22 2 2 4 ... 2 2 ^ n_1 ^G 2 _ 1 G 4 _ 2 ■ ■ • G 2n ) 

< P(4 n _i < - log(e) + (2 n - 1) log(2) + (2 n) log 


Step 2.2. By (1.101 and the arguments developed in Step 2.1, we obtain 
PCM: > n) < P(s n > 2e) < P(Z „_1 - EZ n _ 1 < r/(e,n ) - EZ„-i)> 


where 


r](e, n) := -log(e) + (2 n— l)log(2) + (2n) log y>(0). 

2 V 2 , 


Let us observe that, for any n > 0, m := E[X„] = log(4) + > 0 where p, 

defined by (1.41. Hence 


is 


n k 


E[Z„] = 5>[S„] = = m^(k+l) = 

k =0 k=0 j=0 k —0 


m(n + l)(n + 2) 


Thus, for n large enough, r?(e,n) — EZ n _i < 0. Introducing d n := \mn(n + 
l)/2 — r)(e,n)\, we observe that, for any 0 < k < 1/2 there exists K(k, e) G N 
such that d n > mn 2 { 1/2 — k) for n sufficiently large that is n > K(k, e). After 
straightforward computations, we can choose 


N(«,e) := 


Markov’s inequality leads to 


log(2e)| 


mn 


_ log(2y(0)) 

2k tok 


+ 1 . 


P(M > n) < P(|Z n _i - E[Z„-i]| > d n ) < 


E[(Z n _i - E[M-i]) 4 ] 
dt 


(1.15) 


(1.16) 


Let us note that Xj := Xj — m are i.i.d. random variables with finite moments 

_ 

of any order. We denote mu := E[A ; ], Therefore we obtain 


n— 1 k 


n —1 


IL -L A, ^ 11 i. 

Zn-i :=E[(M-i-E[M-i]) 4 ]=e[(^^M) ] =E[(^(n-i)M) 


k =0 j —0 


3= 0 


n— 1 


= '52(n- j) 4 m 4 + 2 (n- jf(n-kf 

j =0 0<j<k<n —1 


-XK n (n + l)(6n 3 + 9?i 2 + n - 1) + ^ n 2 (n + l) 2 (2n + l) 2 . (1.17) 

ou 00 
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Hence, there exist a constant C 0 > 0 such that E[(Z„_i — E[Z„_i]) 4 ] < Con 6 . 
Combining the previous inequality with (1. 15[) and (1.16) leads to 


P(A/" e > n) < 


Co 


m 4 ( 1/2 — ft) 4 n 2 
Consequently, the following upper-bound holds 

E[A4] = P ( A/ ’ e > n ) ^ £ ) + 


-i, for n>H(ft,e). 


Cn 


E ^ 


n>0 


m 4 (l/2 — ft) 4 z 


In order to conclude, it suffices to note that H(k, e)—>-ooase—5>0, the second 
term in the previous inequality therefore becomes small as e —► 0: the leading 
term is finally H(ft, e) which is equivalent to y/\ log(2e)|/(mft) by (1.15). □ 


2 First-passage time to boundaries with bounded 
derivative 

The algorithm presented in Section |T] is simple to achieve (it only requires inde¬ 
pendent Gaussian random variables) and efficient: the averaged number of steps 
is of the order y/\ log ej where e stands for the small parameter appearing in the 
rejection sampling (see Theorem In order to apply Algorithm |(A1)| the 

curved boundary, the Brownian motion is going to hit, has to satisfies suitable 
conditions: (HI I, (H2) and (|H3|). Asking for the monotonicity of the function p 


is quite restrictive, that’s why we present an extension of the algorithm which 
is of course less efficient (even if the average number of steps is still very small) 
but which permits us to deal with more general boundaries. Let us introduce 
the following assumption: there exist two constants p+ > 0 and p- > 0 such 
that 


ip : K + 


is a (^-continuous function satisfying 


sup < p+ 
t> o 


and 


inf ip' (t ) > — p_. 
t> o 


(H4) 


For such boundaries, we present an algorithm which permits us, for any K £ R + , 


to approximate the hitting time r ^ = t v A K, where t v is defined in (0.11. Let 


us introduce some notations: the inverse Gaussian distribution of parameters 
p, > 0 and A > 0 will be denoted by / (p, A) and is defined by its the probability 
distribution function: 


/(z) = 


A 

2ttx 3 


exp 


[ X(x- p) 2 \ 
l 2p?x i 


L {z>0} • 


Algorithm (A2). Let e > 0 be a small parameter and r > p_ where p- is 
defined in (H4). 


Initialization: (T,H) = (0, <p(0)) and N ti K = 0. 

While H > e and T < K, simulate G an inverse Gaussian random variable with 
distribution I(H/r, H 2 ) and do: 


H i — ip{T + G) — p{T) + r G, 
T t— G + T, 
a f Cj K t— a r e> K +1- 


( 2 . 1 ) 
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Outcome: t* ,k 4— T A K and fif e ,K- 

Algorithm |(A2)| is quite simple, it only requires the simulation of inverse 
Gaussian distributed random variables. Let us recall the following scaling 
property: if G ~ I(H/r, Ft 2 ) then HG/r ~ 1(1, rH). Moreover is 

Chi-squared distributed with one degree of freedom (the square of a standard 
Gaussian random variable). In order to simulate an inverse Gaussian random 
variable, we suggest to use the algorithm introduced by Michael, Schucany and 
Haas (see m or [13 p. 149]). Let us now state the efficiency of Algorithm |(A2)| 
The inverse Gaussian distribution does not permit us to argue in a similar way 
as in Section [T] That’s why we are going to use the general potential theory 
in order to upper-bound the averaged number of steps. This kind of arguments 
was already introduced in convergence results associated to the Random Walk 
on Spheres algorithm which permits the approximation of the solution of the 
Dirichlet problem, see for instance [32]. 


Theorem 2.1. 


1. Let us assume that the boundary function <p satisfies (H4l 


then the random variable t)L ,k defined in Algorithm (A2) 


converges m 


distribution towards r ^ = t v A K where is defined by 123 as e tends 
to zero. More precisely 


F e , K {t-e) - (1+p)\ — < F K (t) < F e K (t), for any t > e, (2.2) 

l/ 7T 


where Fx (resp. F e x) is the cumulative distribution function of rff (resp. 

V ' 

2. There exist positive constants a, b, Kq, Ki and eo such that: for any 
P+ < and any ( K,r ) satisfying (r + kq)K < the random num¬ 
ber of iterations Af e ,x defined in Algorithm (A2) satisfies the following 
upper bound 

E[A4,k] < (a + 6r)|loge|, Ve < e 0 . (2.3) 


3. For non increasing functions ip: there exists two positive constants a and 
eo such that 

E[A4,ic] <ar 2 1<\ loge|, Ve < e 0 . (2.4) 


This theorem is based on the following intermediate statement which is a 
modification of Proposition |1.2| 

Proposition 2.2. Let ( B t , t > 0) be a standard one-dimensional Brownian 
motion. We introduce the following stopping times: Sq = = 0 and for any 

n> 1: 


s n := inf jt > 0 : B t+T K_^ = ~ r ^} ant ^ : = ( s i + • • • + s n ) A K, 

(2.5) 

where the boundary ip satisfies (H4). Then the following properties hold: 

1- (T n K ) n >o is a non-decreasing sequence which almost surely converges to¬ 
wards rff. 
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2. On the event {si + - • - + s n < K}, the probability distribution of s n +\ given 
the a-algebra T n := a{7 } K ,... ,1~ A } is the inverse Gaussian distribution 
I{T~L n /r,T-Ll) with 

'H n ~<p{Tf)-'p{T£ 1 )+rs n . ( 2 . 6 ) 


3. Let M e := inf{n > 1 : <p(T*)-B t k < e}, M K := inf{n > 1, T A = K}, 
Mf = M e A M a . then Tm k and , defined in Algorithm I (A 2 j\ are 


identically distributed, so are Mff andM e ,K- 


Proof of Proposition \2.2\ The first and the third part of the proof are left to 
the reader. They need similar arguments as those presented in Proposition |1.2| 


Here the monotonicity property is just replaced by (H4| which permits us easily 




to prove that T A < . 

Let us now focus our attention to the second part of the statement. Due to the 
definition of s n+1 and since {T A < K}, we get B t k = <p(fTj^-\) ~ rs n- Hence, 
we have 


s n+ i = inf{i > 0 : B t+T K - B t k = p(T*) - B t k - rt} 
= infjf > 0 : W t = 7 i n — rt}, 


where Wt = B t+ j-k — Bj-k is a standard Brownian motion independent of T n 
and the T n adapted r.v. 7 L n is defined by (2.6 1 . The distribution of s„+i cor¬ 


responds to the distribution of the first passage time of the standard Brownian 
motion with drift at the constant level fi n . The probability distribution is well 
known (see, for instance |281 p. 197]): 


P(a„+i G dt\T n ) = / _— exp- 1 


} dt, 


\phrfi ^ 2 1 

we can consequently identify the inverse Gaussian distribution I(fH n /r, hfff). □ 


Proof of Theorem \2. 1\ 

Step 1. We can prove the convergence in distribution of t^ k towards t a using 
similar arguments as those presented in the proof of Theorem |1.1| The upper- 
bound in (2.21 is an adaptation of (1.71 which requires that T A < and that 
and Tm k are identically distributed. These conditions are satisfied, see 


,K 


Proposition 2.2 For the lower-bound in (2.21, we obtain 


F e ,K(t - e) < P(\T m k ~ t a | > e) + F e , K (t), 

r K | 


see (1.8) for the details. Let us note that \t m k — t^\ > e leads to the condition 
tj^k < K. Hence = A4 e . Using the Markov property, the following bound 
holds: 

P(|7^f-<|>e)<l-P( sup B u >e+ sup ip{Tm, + u) - p{Tm, 


v 0<n<e 


0<u<e 


Here ( B t , t > 0) stands for a standard Brownian motion independent of 7 ~M e - 
Combining Hypothesis (H4 ) and the reflection principle of the Brownian motion 
leads to 


P(|7^ e -T*| 


> e) < 1 — P(|H C | > e(l + p .(_)) < (1 + p+) 
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and consequently to the lower bound (2.21. 

Step 2. Let us now focus our attention to the averaged number of steps in 


Algorithm (A2)| denoted by J\f e ,K- A rough description of the method: we aim to 
construct a Markov chain and to describe the associated potential. The classical 
potential theory then permits us to obtain the announced bound. We introduce 
the Markov chain R n := {T rl . %„) for n > 0. We recall that T n = Si + ... + s n is 


defined by (2.5) and hi n by (2.6 1 . The stopping time defined in Proposition 
2.2 can also be interpreted as the first time the Markov chain {R n , n > 0) goes 
out of the domain E := [0, K] x]e, +oo]. 

Let us consider the function f(x,y) = log (y), defined on E, and denote by 
P the infinitesimal generator associated to the Markov chain (R n ) n > o- By 
Proposition 2.2 and for any (t,h) £ E, we obtain 


Pf(t, h) = E log(^(t + G) - <p(t) + r G) 


where G is an inverse Gaussian distributed random variable with the following 
density function: 


p(x) = 


Y2ttx 3 


exp < - 


(h — rxY 


2x 


x > 0. 


By (H4l, ip(t + G) - ip(t) < p+ G, we get 


Pf{t,h) - f{t,h) < log ( 


1 + 




-E 




(2.7) 


Let us find now an explicit upper bound of Pf — /. Using first the change of 
variables u = rx/h and secondly w4 1 /u, we get 


E 



/ o 


>-(x) 


Y 2t 


(h — 

. exp-— dx 

2x 


I hr 
2n , 
I hr 
2n 


log(u) hr(l — u ) 2 , 

exp--- du 


, 3/2 


2 u 


(1 —u)log(u) hr(l — u) 2 , 

exp--- du. 


P / 2 


2 u 


( 2 . 8 ) 


It is then obvious that E log ( 


)»s(¥) 


< 0. Let us now give a more precise upper- 


bound. We set a = hr, then (2.8) emphasizes that E log on ly depends 


on the parameter a, this dependence being continuous. Let us therefore denote 
this function ij>(a) (see Figure [2] below representing Y obtained with the Monte- 
Carlo method sample size: 10 000). 

Simple computations lead to 
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Figure 2: Monte Carlo approximation of the function ip 


V’(a) :=E lo g(-^-)j 


a f ulog(l + u) au 2 , 

— / —- , ' exp — —-- du (2.9) 

2tt ./o (1 + u) 3 / 2 2(1 + w) 


< - 


< - 


< - 


1 


a I' 00 ulog(l + it) au 

2 *1 (! + „)»/> “P-T* 

Z 100 wlog(l + Wa) w , 

V^Jo (a + u>) 3 / 2 ex P-2 dw 

1 Z 100 wlog(l + ru/a) u; , 

V / 27tA/ 2 (a + tu) 3 / 2 2 


Using the inequality (a + u>) < (1 + 2a)w, we get 

log(i + (2a)" 1 ) f x 1 w 

(1 + 2a) 3 / 2 V^ Ji/2 V™ 2 

where G is a standard gaussian r.v. and so P(G > 1/2) w 0.3085 
We deduce from the previous upper-bound that lim a _ >0 + ip (a) = — oo. Moreover 
the right hand side is a non decreasing function with respect to the variable a. 
Hence 

ip{a) < _ log ( 3 / 2 ) P(G > 1/2) m -0.0241, for a < 1. (2.10) 

3v 3 

Let us observe what happens for large values of the variable a. The Laplace 
method implies that 

ip {a) ~-as a —> oo. 

2a 

Let us prove now that there exists a constant c > 0 such that 


V , ( a ) —-, for any a > 1. 

a 


( 2 . 11 ) 
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For a > 1, we get 




< — * 27r 


itlog(l + zt) era 2 

( 1+ „) 3 /2 W - d * 

a f 1 ulog(l + u) au 2 

(1 + u) »/2 exp-—*. 


Due to the convexity property of the logarithm function (log(l + u) > log(2)u) 
and the Cauchy-Schwarz inequality, we obtain 


w “ )£ -^G E[Gf2l “ E[ ° 2l <°>^]) 


< - 


< - 




a2 3 / 2 V2 
log (2) 


a2 3 / 2 


G-GM 


< - 


> ya 

log(2) 

a2 5 / 2 


(1 — \/3e *), for a > 1. 


We deduce that ijj(a) < —c/a with c « 0.0445 when a > 1. Combining both 
inequalities (2.101 and (2.11) leads to the existence of a constant c > 0 such 
that 

1 A (2.12) 


ip[a) < — c A 1^). 


By (2.71, the following upper-bound holds: for f(x,y ) = log(j/), 


Pf(t,h) - f(t, h) < log ( 


1 + 


?±\-r(± 


r / 


(» 


/ P+ 

<-C 




h> 0, t > 0. 


(2.13) 


Due to the definition of p + , we know that 

h < y>(0) V (r + p+)t < <p(0) V (r + p+)K , 


where p is the boundary the process has to hit. In other words, there exist two 
constants kq > 0 and k\ > 0 such that for any p + < kq and any (A', r) satisfying 
(r + Kq)K < K\ the following bound holds | ^ A rj . Hence: 

Pf(t, h) - J(t, h) < A r) =: 

We deduce that the function g(t,h) defined by g(t,h) = 7£(r) ( f(t,h ) — loge) 
satisfies g(t,h ) > 0 for any (t,h) € E and Pg(t,h) — g(t,h) < —1 on E. The 
potential theory therefore implies: 


E[A U,k] < g{ 0, <p(0)) < A.(r)(log((/3(0)) - log(e)). 

We finally deduce the existence of a > 0 and b > 0 such that E[A4,i<r] < 
(a + br )| loge| for e small enough. 

For the particular case of a non increasing boundary function it suffices to vanish 
p + in (2.131 and to apply the same arguments of the potential theory in order 
to get (|2.4 ) □ 
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3 Examples and numerics. 

In this section, we present three different examples which nicely illustrate the 
efficiency of these new algorithms |(A1)| and |(A2)| 


3.1 Brownian hitting time of ip(t) = y/l + cat 

Let us first consider an application of Theorem PL We observe that ip(t) = 
vd + at is an increasing function satisfying (HI), (H2) and (H3l for a £ [0,1]. 
Consequently Algorithm |(Al)| converges and permits us to obtain an approxi¬ 
mation of the hitting time t v . In the figures, we present the link between the 
averaged number of steps and e which characterizes the approximation error 
size. 

The first figure (resp. the second one) concerns: a = 1 (resp. a = 0.01), 
e = 0.5" (n is represented on the horizontal axis) and the number of simulation 
in order to estimate the averaged number of steps is 10 000. 




(a) a = 1 (b) a = 0.01 

Figure 3: E(A4): mean number of steps for e = 0.5" as a function of n. The 
boundary is ip(t) = \/l + at. 

Let us now present the approximate distribution of the hitting time. 




(a) a = 1 


(b) a = 0.01 


Figure 4: Empirical distribution of the approximate first hitting time of the 
boundary tp(t) = s/l + at. 
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3.2 Brownian hitting time of cp(t ) — ct + (3 cos (iot) 


Let us now consider the first time the Brownian motion hits the periodic bound¬ 
ary ip(t) = a -f /3cos(wt). Since the boundary is not an increasing function, we 
shall use Algorithm |(A2)| Theorem 2.1 ensures that the algorithm converges. 
Let us therefore use the Monte-Carlo method in order to estimate precisely the 
average number of steps. As explained in the previous section, the simulation 
procedure permits to approximation of the stopping time !\K for some given 
fixed time K. Figure |5] illustrates the approximation t v by r^ ,A , where the 
parameters are fixed at a = 3.5, (3 = 3 and to = tt/2. The maximal time are 
K = 20 on one hand and I\ = 100 on the other hand and the error rate is given 
by e = 0.5", for 1 < n < 10. A sample of 10A8 paths has been simulated to 
approximate the mean. We know that the mean number of steps is a decreasing 




n t 

(a) E(Wi /2 ",k) versus n (b) Distribution of n = 10, K = 20. 


Figure 5: Approximation of t v with tp(t) = 3.5 + 3cos(7rf/2) 

function of e and an increasing function of K. Figure [6] gives the evolution of 
the mean number of steps as a function of the truncation K. In practice, we 



Figure 6: E(A ^’ K ) as a function of K. 


obtained easily an impressively accurate approximation of t v . 
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3.3 The first time the Ornstein Uhlenbeck process hits the 
boundary ip(t) = a + /3 cos (ut) 

This last section concerns a particular framework where passage times play a 
crucial role: the spiking neuron analysis. Let us roughly explain how neuronal 
firing activities has been modeled. The potential difference that exists across 
the cell membrane is modeled as an Ornstein-Uhlenbeck process. As soon as 
this membrane potential exceeds a given threshold, the neuron releases a rapid 
electrical signal called a spike and the membrane potential is directly reset to 
an initial voltage. Hence the interspike interval is identify with the first hitting 
time of an OU process whereas the spike train forms a renewal process. Such a 
stochastic leaky integrate and fire (LIF) neuronal model is a good compromise 
between realism and mathematical tractability 

The standard OU process can be adapted when the inputs are time-dependent. 
It especially concerns many situations where the sensory stimuli, like sound, con¬ 
tain an oscillatory component. We observe then oscillating membrane potentials 
in the neuron, generating rhythmic spiking patterns (see Iolov, Ditlevsen and 
Longtin m and references therein). 

In such a model with time-dependent forcing, the membrane voltage denoted 
by (V t , t > 0 ) satisfies the following stochastic differential equation: 

dV t = (x{t) - ^ 7 ) dt + a dB t , 


until it reaches the voltage threshold U t h- Here \ represents a current acting on 
the cell, r is the membrane time constant, a is the strength of the stochastic 
fluctuations and finally ( B t , t > 0) stands for the standard Brownian motion. 
The length of the interspike interval is therefore directly related to the passage 
time of the stochastic process (Vt) through the threshold Vth- 

Let us introduce a simple change of variable, given by X t = Vt — ip(t) with 
ip the deterministic function satisfying: 


v'{t) = x(t) 


<f{t) 


¥>( 0 ) 


Vo. 


By straightforward computations, we can prove that X t is a classical OU pro¬ 
cess. In other words, the first passage time of the voltage V t through the given 
threshold V t h is almost surely equal to the first passage time of the OU pro¬ 
cess ( Xt ) through the curved boundary Simulating hitting times to curved 
boundaries for OU processes is therefore a main task. 

That’s why, we focus our attention to the last example which concerns the 
one-dimensional Ornstein-Uhlenbeck process defined by: 


dXt — dBt — A Xt dt 7 Aq — Xq ■ 


(3.1) 


The aim is to approximate the first passage time through the particular simple 
curved boundary ip(t) = a + (3 cos (cot) where <£>( 0 ) > Xq. 

Since the Ornstein-Uhlenbeck process can be represented as a time-changed 
Brownian motion, the question is directly related to the main results of this 


study. Indeed the solution of (3.1) is given by 


X t =e 


-At 


X 0 + 


_A s 


dB. 


1 v — 


t > 0 . 
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Using Levy’s theorem, (X t , t > 0) has the same distribution as (Y t , t > 0) 
defined by 

Y t := e~ xt (x 0 + W„ (t) ) , t > 0, 

with u(t) := -^ ( e 2Xt — 1) and W a standard Brownian motion. We deduce that 

Tp := inf{f > 0 : X t = ip{t)} 
has the same distribution as 

% := inf jt > 0 : e~ xt (x 0 + W u = <£>(t)j 

= inf > 0 : W s = tp^r 1 {s))e Xu ~^ s) - x 0 } 


where 


:= inf{t > 0 : W t = z/(i)}, := \J\ A 2A t <p(^ 


log(l + 2A t) 
2A 


x 0 . 


Consequently, in order to simulate the Ornstein-Uhlenbeck hitting time Tp A I\ 
for some K, we simply use Algorithm |(A2)| and propose an approximation of 
the Brownian hitting time A K with K := u(K) = (e 2XK — 1)/(2A). 

Let us note that a straightforward computation leads to the following upper- 
bound: 

. ,. .. Act: A A 3 A w/3 ^ „ 

\i> {t)\ < - 7= —< Aa A A/3 A lu/3, t > 0. 

V 1 + ZAt 


In other words, the continuous curve i^> satisfies Hypothesis (H4l: Algorithm 
|(A2)| therefore converges and Theorem 2.1 can be applied. 

In the following numerical experiences, we will choose r = 0.5 A Act: A A/3 A w/3. 
Figures [7] and [8] concern the following choice of parameters: Xo = 0, a = 2, 
/3 = 1, w = 7 t/ 5, A = 0.5. We have chosen K = 5 for Figure [7] and K = 10 for 
Figure [8j In both cases, the first figure represents the average number of steps 
as a function of n where the approximation parameter e is chosen as 0.5™, for 
n = 1, • • • ,10. The average has been estimated using 5.10A6 simulations. The 
second figure represents the distribution of Tp A K for n = 10. 

We observe that the change of time K = (e 2AA — 1)/(2A) increases very fast 
with K and the number becomes quite large when K increases. Note however 
that the number of random variables we have to simulate keeps relatively small 
in comparaison with the use of a classical stopped Euler scheme usually used 
to approximate Tp. Each Euler scheme introduces a bias, which goes to 0 as 
the time discretization length goes to 0. We have plotted the error on Figure [9] 
for the approximation of K(T^). We can observe that a simple Euler scheme 
has an order of convergence 1/2. We illustrate the improvement of this rate 
of convergence using two particular modifications of the scheme: the first one 
from [7-fl (i.e. we take into account the first order term of the probability to 
hit the boundary between two time successive time-steps), the second-one from 
|25| (i.e. an adapted modification of the boundary). Both modifications yield 
to a scheme of order 1. 

The numerical cost of our algorithm increases very slowly as the parameter 
e goes to 0. The numerical comparisons are done with e = 2 -20 , such that the 
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error is almost negligible. The time we need is similar to the time for an Euler 
scheme with step 0.01 and the Brownian bridge modification with time step 
0.02. Empirically, we conclude that our scheme over performs previous ones if 
one needs an accuracy larger than those obtained with an Euler scheme with 
time step 0.01. 




(a) E(.V'| / 2 ‘" : k ) versus n 


(b) Distribution of T^’ <L (e = 1/2 10 ). 


Figure 7: First hitting time of ip(t) = a + (3 cos(cot) by an Ornstein Uhlenbeck 
process solution of (|3.1|) (a = 2, /3 = 1, to = 2i r, A = 0.5, K = 5.) 




(a) E(A/*i/ 2 n ,K) versus n 


(b) Distribution of (e = 1/2 10 ). 


Figure 8: First hitting time of ip{t) = a + /? cos (cut) by an Ornstein Uhlenbeck 
process solution of ( |3.1| ) (a = 2, = 1, uj = 2i r, A = 0.5, K = 10). 
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Figure 9: Bias for Euler schemes to approximate E(T^), the mean first hitting 
time of <p(t) = a + p cos (uit) by an Ornstein Uhlenbeck process solution of (3.1) 
(a = 2, P = 1, w = 2tt, A = 0.5, K = 5). 
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